options(digits = 3)
p <-
  brazil_data %>% 
  filter(urb == 1) %>%
  select(pop_2010, private_center_prop, private_school_prop) %>%
  mutate(private_school_prop = private_school_prop*100) %>%
  pivot_longer(cols = private_center_prop:private_school_prop, names_to = "names", values_to = "values") %>%
  mutate(condition_label = case_when(names == "private_center_prop"  ~ "Private Health Centers (%)",
                                     names == "private_school_prop" ~ "Private School Enrollment (%)")) %>%
  ggplot(aes(x=log10(pop_2010), y = values)) + 
  geom_point(color="darkgrey", shape = 21, alpha = 0.5) + 
  geom_smooth(method = "loess", color = "black") +
  facet_wrap(condition_label ~., scales = "free") +
  stat_summary_bin(fun='mean', bins=20,
                   shape = 21, fill = "lightgrey",
                   color='black', size=2.5, geom='point') +
  theme_bw() +
  scale_colour_grey() +
  scale_x_continuous(limits = c(3, 7), 
                     breaks = c(3, 4, 5, 6, 7),
                     labels = c("1,000", "10,000", "100,000", "1,000,000", "10,000,000")) +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.line.y.left = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        legend.title = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.title.y = element_blank()) +
  xlab("Population, log scale")




ggsave("./_4_outputs/figure_7.pdf", plot = p, width = 10, height = 4)
